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Abstract 


We  present  a  method  to  determine  whether  a  set  of  equations  has  a  non-negative  integer  solution.  The  method  is  designed 
for  the  particular  occurence  of  this  problem  in  the  context  of  compiler  analysis  of  parallel  programs.  The  system  of 
equations  is  first  transformed  to  Smith  normal  form  to  determine  if  any  integer  solutions  exist.  In  case  of  multiple  integer 
solutions,  a  parameterized  solution  space  representing  all  non-negative  solutions  is  obtained.  Fourier-Motzkin  elimination 
is  employed  to  determine  if  the  real  solution  space  is  empty.  If  the  solution  space  is  not  empty,  either  the  existence  of  an 
integer  solution  is  readily  verified,  or  a  simplified  convex  region  is  obtained  such  that  the  original  system  of  equations 
has  a  solution  if  and  only  if  this  convex  region  contains  an  integer  point.  The  main  result  of  the  paper  is  a  set  of  new 
heuristic  search  procedures  that  are  used  to  identify  an  integer  solution  in  a  convex  region,  or  to  prove  that  no  integer 
solution  exists.  These  are  based  on  the  geometrical  properties  of  convex  regions  that  are  not  empty  but  do  not  contain 
integer  points.  This  method  is  an  exact  and  efficient  way  of  solving  integer  programming  problems  arising  in  dependence 
and  synchronization  analysis  of  parallel  programs. 
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1  Introduction 


Sever ai  mathematical  problems  that  arise  in  the  development  of  compilers  for  parallel  languages  can  be  transformed  to 
Integer  programming  problems.  Determining  whether  a  data  dependence  exists  between  two  array  references  can  be 
transformed  to  the  problem  of  determining  whether  an  Integer  solution  to  a  set  of  equalities  and  inequalities  exists  (1,10]. 
Similarly,  in  an  event  variable  synchronization  model,  determining  whether  the  synchronization  present  in  a  program  is 
sufficient  to  protect  a  dependence  can  be  transformed  to  the  problem  of  determining  whether  a  system  of  linear  equations 
has  a  non-negative  Integer  solution  [2, 9].  Both  these  problems  require  integer  programming  which  is  known  to  be  NP- 
complete.  In  this  report,  we  present  a  practical  solution  technique  for  determining  whether  a  non-negative  integer  solution 
for  a  system  of  equations  exists.  This  is  equivalent  to  the  transformed  synchronization  analysis  problem  mentioned  above, 
and  the  solution  methodology  can  also  be  used  to  solve  the  data  dependence  analysis  problem.  Both  these  problems 
have  characteristics  that  make  it  possible  to  develop  methods  that  work  well  in  practice,  even  though  the  problems  are 
NP-complete.  We  state  the  characteristics  of  systems  of  equations  that  are  constructed  from  the  synchronization  analysis 
problem: 

•  The  number  of  variables  corresponds  to  the  number  of  synchronization  operations  relevant  to  a  specific  data 
dependence,  and  is  expected  to  be  small. 

•  The  number  of  equations  corresponds  to  the  number  of  subscript  positions  in  the  array  references  causing  a  data 
dependence,  and  is  expected  to  be  small. 

•  The  coefficients  of  terms  are  synchronization  distance  vector  components  and  are  usually  small  integers,  very  often 
one  or  zero. 

Systems  of  equations  constructed  from  the  dependence  analysis  proble  have  similiar  characteristics.  The  point  is  that 
these  equation  systems  are  fairly  small  and  can  be  solved  efficiently  in  practice. 

We  outline  the  main  steps  in  the  solution  method.  We  first  solve  the  equation  system  for  all  (not  necessarily  non¬ 
negative)  integer  solutions.  This  is  often  referred  to  as  solving  linear  diophantine  equations.  If  such  a  solution  is  unique 
or  does  not  exist,  we  determine  this  fact  and  do  not  need  to  proceed  any  further.  If  there  are  multiple  solutions,  we  obtain 
parametric  equations  that  describe  the  solution  space.  Next  we  determine  the  non-negative  solution  space  and  check 
whether  it  is  empty.  If  it  is  not  empty,  we  use  heuristic  search  procedures  to  determine  if  it  contains  an  integer  solution. 


2  Solving  Systems  of  Linear  Diophantine  Equations 

The  theory  of  solving  linear  diophantine  equations  has  been  discussed  in  many  texts  addressing  number  theory  and  integer 
programming  [5, 4],  Our  contribution  is  to  present  a  solution  method  in  a  way  that  can  be  easily  programmed  on  a 
computer  and  analyze  the  complexity.  We  also  discuss  why  we  chose  the  specific  method  that  we  present. 

The  problem  can  be  stated  as  follows: 

Find  all  integer  solutions  of 

Mlmxn'M"*1  =  Wmxl 

We  first  present  the  necessary  mathematical  concepts,  and  then  an  algorithm  to  solve  the  problem. 

2.1  Mathematical  Concepts 

Following  are  the  elementary  column  operations  on  a  matrix: 

1.  exchanging  two  columns 

2.  multiplying  a  column  by -1 

3.  adding  an  integral  multiple  of  one  column  to  another  column 

Corresponding  elementary  row  operations  are  similarly  defined.  We  shall  use  elementary  operations  to  transform  a  matrix 
to  a  form  desirable  for  solving  the  equations.  In  the  rest  of  this  section,  we  introduce  two  kinds  of  matrices  that  describe 
specific  elementary  operation  sequences  that  we  shall  use. 

A  transposition  matrix  [P]„xn  is  a  square  matrix  of  Is  and  Os  in  which: 
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1.  each  row  and  each  column  has  exactly  one  1. 

2.  all  the  l's  except  two  are  along  the  main  diagonal. 

The  matrix  shown  below,  Pu  is  a  transposition  matrix  of  order  5.  The  subscript  25  signifies  that  the  1$  in  rows(and 
columns)  2  and  5  are  out  of  place  relative  to  a  unit  matrix. 


Pu  = 


1  0  0  0  0' 

0  0  0  0  1 

0  0  10  0 

0  0  0  1  0 

0  10  0  0 


Theorem  1  Premultiplying  a  matrix  A  with  a  transposition  matrix  Pa  is  equivalent  to  interchanging  rows  i  and  j  of  A, 
and  postmultiplying  A  with  Py  is  equivalent  to  interchanging  columns  i  and  j  of  A. 


We  show  this  with  an  example  where  premuitiplying  A  with  Pis  yields  A  with  rows  2  and  5  permuted. 
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A  Subtraction  matrix  is  a  square  matrix  in  which: 

1.  all  items  along  the  main  diagonal  are  1. 

2.  S[t,  j]  =  —a. 

3.  all  other  elements  are  0. 

Following  is  an  order  5  subtraction  matrix  Si3,a 


$13, a 


1  0  -a  0  0 

0  10  0  0 
0  0  10  0 
0  0  0  1  0 

0  0  0  0  1 


Theorem  2  Premultiplying  a  matrix  A  with  a  subtraction  matrix  Sijia  has  the  effect  of  reducing  the  ith  row  of  A  by  a 
times  the  jth  row  of  A.  Postmultiplying  a  matrix  A  with  a  subtraction  matrix  5^,0  has  the  effect  of  reducing  the  jth 
column  of  A  by  a  times  the  ith  column  of  A 


We  show  this  with  an  example  where  postmultiplying  A  with  5 1  j,0  yields  A  with  every  element  of  column  1  reduced 
by  a  times  the  corresponding  row  element  in  column  3. 
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Multiplication  of  two  subtraction  matrices,  say  Siy,01  and  Sn,a ,  yields  a  matrix  with  la  along  the  diagonal,  <*i  and  a2 
at  locations  Sy  and  S,*  respectively,  and  0s  elsewhere.  The  operation  of  subtracting  different  multiples  of  a  column  of 
a  matrix  from  other  columns,  which  can  be  accomplished  by  successive  multiplication  by  different  subtraction  matrices, 
can  be  achieved  by  a  multiplication  with  a  composite  subtraction  matrix  constructed  as  described  above. 

Transposition  matrices  and  subtraction  matrices  are  examples  of  a  more  general  class  of  matrices  called  regular 
unimodular  matrices.  A  regular  unimodular  matrix  is  a  square  matrix  whose  determinant  has  the  value  (+ 1 )  or  ( - 1 ) .  We 
state  some  of  the  properties  of  regular  unimodular  matrices  that  are  relevant  for  our  purposes: 
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•  the  product  of  two  regular  unimodular  matrices  Is  a  regular  unimodular  matrix. 

•  the  inverse  of  a  regular  unimodular  matrix  formed  of  all  integers  always  exists  and  is  a  regular  unimodular  matrix. 


2*2  Smith  Normal  Form 


In  this  section  we  define  a  matrix  form  called  the  Smith  normal  form  [8]  and  how  a  matrix  can  be  transformed  to  Smith 
norma}  form  using  operations  defined  by  transposition  and  subtraction  matrices.  The  theory  of  transforming  matrices 
to  Smith  normal  form  and  solving  systems  of  equations  in  integers  using  it  can  be  found  in  [5, 4}.  Here  we  present  an 
algorithmic  description  of  the  method  and  argue  that  it  is  a  suitable  first  step  for  solving  the  systems  of  equations  we 
expea  to  encounter. 

The  Smith  normal  form  of  a  matrix  [j4]mXl»  of  rank  r  is  a  matrix  D  of  the  following  form: 


di 

0 

...  o 

0  ••• 

0  ' 
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0  ... 

0 
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•  • 

...  dr 

0  ... 
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...  o 

0  ... 

0 

.  0 

0 

...  o 

0  ... 

0  . 

where  dt  ^0,  k  =  1, 2,  •  •  •  r. 

Theorem  3  For  every  matrix  [A]mxn,  there  exist  regular  unimodular  matrices  [f/]mxm  aru/[V]nxn  such  that 

[U]  mxm  •  M  mxn*  [V]  nxn  ~  [D\  mxn 


where  [D]mXn  is  in  Smith  normal  form 

We  present  a  procedure  to  compute  the  Smith  normal  form  of  a  matrix  and  the  corresponding  unimodular  multiplier 
matrices. 

algorithm 

Input:  Matrix  [X]mx„ 

Output:  Matrices  [D]mxn,  [U]mxm  and  [V]„xn  such  that 

mxn*  [V)nxn  =  [D]  mxn 


where  U  and  V  are  regular  unimodular  matrices  and  D  is  in  Smith  normal  form. 

Initialization:  D  =  A,  U  and  V  are  unit  square  matrices  of  dimensions  m  and  n  respectively. 

Method :  We  use  the  procedure  FindSmith  stated  in  Figure  1,  which  converts  a  matrix  of  the  form  Dr  to  Dr+, 
below  in  each  step. 


di  0 
0  ^2 
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[0]r 


xn — r 
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0  0  ...  dr 


[0]m-rxr  [J^jm-rxn-r  „ 
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[0]m-'r-lxr+!  (^4. l]m*-r-!xn-r-l  J 

The  procedure  stops  when  0  is  transformed  to  a  Smith  normal  form. 
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procedure  FindSraith(D) 


r  *  1 

U  and  V  are  unit  matrices  of  rank  m  and  n  respectively. 

aize(D)  is  the  minimum  of  numbers  of  rows  and  number  of  columns  in  D. 

We  refer  to  the  submatrix  of  D  formed  of  elements  with  row  and  column  number 
greater  than  r  as  D£ . 

1.  If  r  >  aize(D)  or  the  matrix  UT  has  all  Os  STOP. 

/*  D  is  already  in  Smith  normal  form  V 

2.  Determine 

minnow,  mincol  and  minval  -  the  location  and  value  of  the  smallest  nonzero  number 
in  D1,.. 

3.  Interchange  row  r  with  minnow  and  column  r  with  mincol. 
f*  Now  location  Drr  has  value  minval.  */ 

Premultiply  V  and  postmultiply  U  with  permutation  matriens  that  reflect 
this  transformation. 

V  *  Pminrow ,r  'V  and  U  ~  U. Pmincol ,r 

4.  Reduce  row  and  column  r  with  appropriate  subtraction  matrices: 

For  s'  ss  r  +  l,n: 

ai  =  Dri  div  minval 
Dri  =  Dri  mod  minval 

The  relevant  composite  subtraction 

matrix  5  is  an  n  x  n  unit  matrix  except  that  ar+1, arr+i, •  •  •  a„  are  the  entries  in 
row  r,  columns  r+l  through  m,  respectively.  (see  earlier  discussion). 

V  =  S.V 

Repeat  the  subtraction  procedure  for  column  r  and  postmultiply  U 
accordingly. 

5.  If  row  and  column  r  of  D  contain  all  Os,  r  =  r  +  1  GOTO  step  2.  Else  GOTO 
step  1. 

end  of  procedure 


Figure  1:  Conversion  of  a  Matrix  to  Smith  normal  form 


Complexity:  We  analyze  the  complexity  of  procedure  FindSmith  shown  in  Figure  1.  Steps  1  and  2  search  a  submatrix 
of  [Ojmxn  in  O(mn)  time.  Step  3  interchanges  one  pair  of  rows  and  one  pair  of  columns  and  multiplies  corresponding 
permutation  matrices  to  unimodular  matrices  U  and  V,  which  is  again  equivalent  to  a  row  pair  and  a  column  pair 
interchange.  Hie  time  taken  is  0(m  +  n).  Step  4  performs  addition  of  a  multiple  of  a  column  to  the  remaining  columns 
in  a  submatrix  of  D  and  a  similar  operation  with  rows  of  D.  A  similar  operation  is  performed  on  unimodular  matrices  U 
and  V,  when  they  are  multiplied  by  the  subtraction  matrix  computed  above.  Ail  of  these  operations  take  O(mn)  time. 
Thus  the  overall  complexity  of  one  execution  of  steps  1-4  is  0(mn).  • 

In  step  5,  the  control  goes  back  to  step  2  if  row  and  column  r  do  not  have  all  Os,  except  in  location  Drr.  If  Drr 
became  1,  the  above  condition  must  be  met  in  the  next  step.  Also  Drr  is  reduced  in  every  step,  in  a  manner  similar  to 
computing  ged  of  a  set  of  numbers  and  is  expected  to  reach  L  in  log  t  (value)  steps,  where  value  is  the  smallest  element  in 
the  submatrix  D'r  at  the  beginning  of  this  reduction  step.  Finally,  the  loop  from  1-5  is  executed  R  times,  where  R  is  the 
rank  of  the  original  matrix  A. 

Thus  the  complexity  of  computing  the  Smith  normal  form  and  the  corresponding  unimodular  matrices  for  matrix 
[■AJmxn  of  rank  R  is  0(mn/Hog(va/tie))  where  value  can  be  approximated  by  the  median  of  non-zero  values  in  matrix 
A.  Thus  if  aize  is  the  higher  dimension  of  A,  the  complexity  of  the  operation  is  bounded  by  0((«re)3log(t>aiue)). 
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cad  of  algorithm 

Systems  of  equations  in  integers  can  also  be  solved  by  transforming  a  matrix  to  a  triangular  matrix  form  in  which  only 
row  or  column  operations  are  used  in  the  transformation  procedure.  However,  in  those  methods  a  new  lowest  element 
in  pivot  position([£7]rr  in  earlier  discussion)  can  be  chosen  only  from  a  particular  row  or  column.  In  computing  Smith 
normal  form,  the  new  pivot  can  be  chosen  from  anywhere  in  the  unprocessed  section  of  the  matrix,  hence  there  is  a  greater 
chance  of  finding  a  pivot  that  would  zero  out  the  corresponding  row  and  column. 


23  Solving  Equation  Systems 


We  show  how  transformation  to  Smith  normal  form  is  used  to  solve  a  system  of  equations  in  integers. 

Let 

M]mxn-[*]nxl  —  Wmxl  (1) 


be  a  system  of  linear  equations  for  which  we  need  to  determine  the  integer  solution  set. 

Let 

[®!mxn  = 

(2) 

where  D  is  in  Smith  normal  form  and  U  and  V  are  unimodular  matrices. 

We  rewrite  from  equation  1 

[y]mXm-[-<4]mxn-[V]nxn'[V  *]nxn-[*]nxl  =  [V]mxm  ^]mxl 

(3) 

l^]mxn't^lnxn'[*]nxl  =  [I7]mxm[&]mxl 

(4) 

We  define 

[f]»xi  =([V-‘].[*])„X, 

(5) 

From  equation  4  we  get 

[I7]mxn  •[*]»*  xl  =‘[^]mxm[&]mxl 

(6) 

If  r  is  the  rank  of  matrix  A,  and  hence  of  matrix  D,  then  [DmXn]  is  of  the  form  f  ^r*T  j? 

.  We  rewrite  equation  6 

as  two  equations  representing  the  top  r  rows  and  bottom  m  -  r  rows  of  the  column  matrices  being  equated. 

[D;Ur.[trrUl  =  Kim-Mfflrxl  • 

(7) 

[0](m-r)x,  =  t«y][6!)^:rr](m-r)x, 

(8) 

Thus,  for  the  original  system  of  equations  to  have  an  integer  solution 

KJr*l  =  (Orr]r-‘r(([f/].[6]);]rxl 

(9) 

must  be  formed  of  integers  and 

K(^M6])mIr](">-r)xl  =  [0j(m-r)xl 

(10) 

If  an  integer  solution  exists,  it  can  be  represented  by  rewriting  from  equation  4. 

[*jnxl  =  (I^Jn  xi»'[f]n  X 1 

(11) 

We  demonstrate  the  complete  procedure  of  solving  systems  of  equations  in  integers  with  an  example.  Consider  die 
following  system  of  equations: 


We  restate  the  system  as: 


2xi  +  3*2  +  4*3  +  5*4  —  3 

*i  -  *2  +  *3  +  2*4  =  5 

2*i  +  0*2  +  2*3  +  5*4  =  -3 


'  2  3  4  5 

1-112 
2  0  2  5 


*i 

*2 

*3 

*4 
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We  have 


‘2  3  4  5' 

'10  0' 

10  0  0' 
0  10  0 

0  0  10 

II 

o 

Q 

II 

1-112 

2  0  2  5 

U°  = 

0  1  0 

0  0  1 

0  0  0  1 

We  shall  transform  these  matrices  such  that 

Dk  = U".A.Vi 

remains  an  invariant 

First  we  interchange  rows  I  and  2  of  D  to  obtain  the  smallest  non-zero  element  of  the  matrix  in  the  top  left  corner.  In  the 
above  equation  we  premultiply  D  and  U  by  permutation  matrix  Pit-  We  have 


'  1  -1  1  2  " 

'010' 

Dl  = 

2  3  4  5 

f/'  = 

1  0  0 

2  0  2  5 

0  0  1 

To  get  the  smallest  possible  values  in  1st  column  of  D,  we  premultiply  by  composite  subtraction  matrix  Sti.t- 
For  the  same  effect  on  first  row,  we  postmultiply  by  .Sj3,t  .Su^  We  have 


’  1 

0 

0  ' 

SpRB  — 

-2 

1 

0 

■Spost  - 

-2 

0 

1 

After  corresponding  multiplications  we  get 


1  1  -1 

0  i  0 

0  0  1 

0  0  0 


-2  ' 
0 
0 
l 


1  0  0  0" 

'  0  10' 

D2  = 

0  5  2  1 

0  2  0  1 

u3  = 

i  -2  0 

0  -2  1 

1  -1 

1  0 

0  1 

0  0 


-2  1 


0 

0 

1 


Again  to  move  the  lowest  nonzero  element  of  D  to  D#  we  interchange  columns  2  and  4  of  D.  We  postmultiply  D  by 
permutation  matrix  Pu.  We  get 


'  1  0  0  0' 

'  1  -2  -1  1  ' 

DJ  - 

0  12  5 

0  10  2 

f/3  =  V2  v>  = 

0  0  0  1 

0  0  10 

0  10  0 

To  reduce  the  row  and  column  2  of  matrix  D,  we  premultiply  by  subtrarion  matrix  Sn.t  and  postmultiply  by  composite 
subtraction  matrix  Sts, i-Sm.j  and  obtain 


’  1 

0 

0 

0  ' 

0 

1 

0  ' 

'  I 

-2 

3 

11  ' 

D*  = 

0 

1 

0 

0 

u*  = 

1 

-2 

0 

V*  = 

0 

0 

0 

1 

0 

0 

-2 

-3 

-1 

0 

1 

0 

0 

1 

0 

L 

0 

1 

-2 

-5 

Now  the  lowest  nonzero  element  in  the  part  of  D  not  in  Smith  normal  form  is  already  in  the  desired  position.  To 
reduce  the  remaining  one  element,  we  postmultiply  by  subtraction  matrix  Ss*,  i  and  get: 


'  1 

0 

0 

0  ' 

'  1 

-2 

3 

8  ' 

D5  = 

0 

1 

0 

0 

Vs  =  U*  Vs  = 

0 

0 

0 

1 

0 

0 

-2 

-l 

0 

0 

1 

-1 

0 

1 

-2 

-3 
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We  interchange  columns  3  and  4  to  get  the  smallest  element  in  the  desired  position  by  postmultiplyingby  permutation 
matrix  P54  and  gee 

1  0  0  0  ‘ 

0  1  0  0  U6  =  US  V6  = 

0  0-1  -2 

Postmultiplyingby  subtraction  matrix  Su.i  brings  D  in  Smith  normal  form  and  we  get  the  final  forms: 


‘  1  0  0  0  ]  [0  10' 

D=  0  1  0  0  U  =  1-2  0  1/  = 

0  0  -1  0  j  -1  0  1 


And  we  verify  that 


D  —  U.A.V 


The  top  left  square  submatrix  of  D  with  number  of  rows  equal  to  rank  r  of  D  is  the  following  matrix: 


[D, rJjxJ  = 


10  0  ‘ 
0  I  0 
0  0-1 


Its  inverse  can  be  computed  easily  since  it  is  a  diagonal  matrix. 

[10  0  ' 
=01  0 
0  0  -1 

From  equation  9  we  have 


Since  it  evaluates  to  an  integer  column  and  equation  10  is  trivially  satisfied,  the  system  must  have  integer  solutions.  The 
solutions  of  the  system  can  be  represented  using  equation  1 1 


67  -  13<« 
6  -  2f« 
-6  +  3<4 
-25  +  4t4 


We  rearrange  the  solution  and  replace  (4  by  I  since  there  is  only  one  parameter. 


*i  =  -  I3f  +  67 

=  -2t  +  6 

ij  =  3*  —  6 

*4  =  4(  —  25 


2.4  Comparison  with  Generalized  GCD  test 

Generalized  GCD  Test  is  a  well  known  way  of  transforming  a  system  of  equality  constraints  to  a  system  of  inequality 
constraints.  Both  methods  essentially  use  matrix  elimination  procedures  in  the  domain  of  integer  solutions.  In  our  method, 
any  matrix  element  can  be  chosen  as  pivot,  and  one  row/column  is  eliminated  when  the  pivot  element  is  a  factor  of  every 
element  in  the  corresponding  row  and  column  (trivial  if  the  absolute  value  of  pivot  is  1 ).  In  generalized  GCD  Test,  only 
elementary  row  operations  are  used,  hence  the  choice  of  a  pivot  is  restricted  to  a  single  column.  Thus,  although  we  spend 
more  time  looking  for  a  pivot  element,  we  expect  perform  less  steps  to  eliminate  a  row/column. 
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3  Finding  Solutions  in  a  Feasible  Region 

We  first  summarize  the  results  of  the  solution  procedure  of  the  last  section,  in  the  context  of  the  overall  problem.  We 
started  with  a  system  of  m  equations  (and  say  r  independent  equations)  in  n  variables.  If  r  is  greater  than  or  equal  to  n, 
any  integer  solution  that  exists  must  be  unique,  and  we  would  have  already  determined  whether  a  non-negative  integer 
solution  exists.  However,  if  r  Is  less  than  n,  then  we  have  all  possible  solutions  for  the  n  variables  parameterized  by 
n  -  r  new  variables  which  can  take  any  integer  value.  We  need  to  determine  whether  there  are  any  non-negative  integer 
solutions  to  the  original  problem.  For  the  example  in  the  last  section,  the  problem  is  parameterized  in  terms  of  the  new 
variable  t,  and  the  original  problem  has  a  non-negative  integer  solution  if  and  only  If  the  system: 

— 13<  +  67  >  0 

-21  +  6  >  0 

3  i  -  6  >  0 

41  -  25  >  0 

is  feasible  for  some  integer  value  of  1. 

This  inequality  system  is  particularly  easy  to  solve  since  there  is  only  one  variable  involved.  The  above  system 
simplifies  to: 


t 

< 

67/13 

(1) 

t 

< 

3 

Q) 

t 

> 

2 

(3) 

t 

> 

25/4 

(41 

The  system  is  infeasible  since  (2)  and  (4)  cannot  be  simultaneously  satisfied. 

3.1  Fourier-Motzkin  Elimination  Method 

The  number  of  variables  in  a  system  of  inequalities  can  be  successively  reduced  by  using  Fourier-Motzkin  elimination 
method.  Any  system  of  inequalities  can  then  be  solved  by  successive  reduction  and  back  substitution.  In  this  method,  a 
variable  is  eliminated  by  creating  a  new  inequality  for  each  pair  of  inequalities  in  which  the  variable  being  eliminated  has 
a  different  sign.  This  is  an  important  step  in  the  solution  procedure,  but  our  approach  is  as  same  as  that  taken  in  several 
other  systems.  However,  for  completeness,  wee  will  illustrate  this  method  with  an  example.  A  complete  treatment  can 
be  found  io  [3].  Consider  the  following  set  of  inequalities: 


Ox  -  y  +  6 

> 

0 

(1) 

x  +  2y-  6 

> 

0 

(2) 

-2x  -  3y  +  7 

> 

0 

(3) 

— 4x  -  5y  -H  40 

> 

0 

(4) 

Suppose  we  decide  to  eliminate  x  from  this  system  of  inequalities1,  x  has  a  positive  sign  in  (2)  above  and  a  negative 
sign  in  (3)  315(1  (4)  above.  By  combining  inequality  pairs  (2)-(3)  and  (2)-(4)  to  eliminate  x,  we  obtain  the  following 
system: 

-y  +  6  >  0  (1) 

y-5  >  0  from  (2)  and  (3)  (2) 

3y  +  16  >  0  from  (2)  and  (4)  (3) 

This  simplifies  to: 

5  <  y  <  6 

1  Heuristic*  for  choosing  «  vsrinbte  so  is  to  minimize  the  work  are  discussed  in  (3|. 
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Thus  a  solution  set  to  the  inequalities  does  exist  We  substitute  the  maximum  or  minimum  possible  value  of  y  back  in  the 
original  equations  so  as  to  obtain  the  range  of  values  z  can  take.  In  the  form  in  which  our  equations  are,  we  substitute 
such  that  the  y  term  has  the  minimum  possible  value.  We  get : 


z  +  4 

> 

0 

(1) 

-2z  -  11 

> 

0 

(2) 

-4  x  +  10 

> 

0 

(3) 

These  are  equivalent  to  the  following  inequality: 

-4  <  x  <  -5/2 

Solving  systems  of  inequalities  by  elimination  tells  u-  whether  a  real  solution  exists  to  the  system  of  inequalities,  in  case 
such  a  solution  does  exist,  we  can  determine  the  range  of  the  values  that  the  variables  may  have  in  the  solution  space.  It 
still  needs  to  be  determined  wbether  an  integer  solution  to  the  system  exists. 

Solving  a  system  of  inequalities  by  Fourier-Motzldn  elimination  techniques  can  potentially  increase  the  number  of 
inequalities  from  n  to  (n/2)J  in  every  stage.  Duffin  [3]  discusses  ways  to  lower  the  computation  effort  needed  to  reach 
a  solution.  They  include  rules  for  selection  of  the  variable  to  be  eliminated  at  each  stage,  and  identification  and  deletion 
of  redundant  inequalities.  He  further  argues  that  with  these  modifications,  Fourier  analysis  is  a  practical  method  to  sc've 
systems  of  inequalities.  The  number  of  variables  and  inequalitites  that  we  expect  in  our  problems  is  small  and  we  expect 
the  analysis  to  be  fast  and  take  only  a  small  number  of  operations. 

4  Identifying  Integer  Solutions 

The  system  of  inequalities  obtained  by  solving  diophantine  equations  defines  a  convex  region  in  n  dimensional  space 
where  n  is  the  member  of  variables  in  the  parametric  equations.  We  are  interested  in  determining  the  feasibility  of 
an  integer  solution,  if  the  Fourier-Motzkin  analysis  shows  that  the  solution  space  is  not  empty.  In  Fourier-Motzkin 
elimination,  we  also  determine  the  bounds  in  each  dimension  for  the  convex  region  defined.  Furthermore  these  bounds 
are  tight  in  the  sense  that  at  least  one  vertex  of  the  convex  region  must  be  on  each  bounding  hyperpiane. 

We  define  some  terms  here  that  we  shall  use  in  the  rest  of  this  section.  An  n  -  rectangle  is  a  closed  _oovex  region 
in  n  dimensions  bounded  by  hyperplanes  that  are  parallel  to  the  axis  hyperplanes  (that  is,  Defined  by  .V  =  constant 
for  some  co-ordinate  axis  X \  An  n-rectangle  tightly,  or  strictly  bounds  a  convex  region  if  the  the  convex  region  is 
completely  contained  in  the  n -rectangle,  and  every  hypetplane  defining  the  n-rectangle  intersects  with  the  convex  region. 
The  centroid  of  an  n-rectangle  is  the  point  with  co-ordinates  equal  to  the  midpoint  of  the  range  of  the  n-rectangle  along 
each  axis.  IfXi,Xi,  ...Xn  are  the  co-ordinate  axes  and  (01,03,  ...a„)and  [f)\,Pu  ../3n)  arc  two  points  inn  space,  a  line 
in  n  dimensions  between  them  is  defined  by: 

-  02  _  Xn  -  On 
Pi  -  £*!  Pl~  02  Pn  -  On 


4.1  The  Case  of  Two  Dimensions 

[f  the  system  of  inequalities  obtained  has  only  one  parametric  variable,  the  system  can  be  trivially  sol  ved  as  shown  earlier. 
We  first  describe  a  solution  procedure  when  there  are  two  parametric  variables.  Subsequently  we  will  describe  the  solution 
method  for  three  or  more  variables  by  solving  a  set  of  two  variable  problems.  We  expect  that  most  practical  situations 
would  yield  a  system  in  two  or  fewer  variables. 

A  system  of  inequalities  in  two  dimensions  describes  a  convex  region  in  two  dimensions.  For  the  present  we  assume 
that  the  region  is  finite.  Fourier-Motzkin  elimination  discussed  in  the  last  section  would  yield  a  tight  bounding  rectangle 
for  the  convex  region.  The  bounding  rectangle  is  known  to  contain  integer  points  and  we  want  to  determine  if  the  convex 
region  itself  has  any  integer  points. 

The  following  theorem  holds  for  convex  regions  in  two  dimensions: 

Theorem  4  If  a  finite  con  vex  region  in  two  dimensions  is  tightly  bounded  by  a  rectangle,  then  the  centroid  of  the  rectangle 
is  always  a  part  of  the  convex  region. 
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Figure  2:  Illustration  for  Theorem  4 


Proof: 

Let  C  be  a  convex  region  in  two  dimensions  and  let  R  be  the  rectangle  that  bounds  it  tightly.  Let  abed  be  the  vertices 
of  the  rectangle  R  as  shown  In  Figure  2.  Let  Peentroid  be  the  centroid  of  R,  which  is  also  the  point  of  intersection  of  the 
diagonals  ac  and  bd.  Since  R  bounds  C  tightly,  each  of  the  sides  ab,  be,  cd  and  da  must  contain  at  least  one  point  that 
belongs  to  C.  Let  Ptt*  and  P%c  be  points  on  the  sides  ab  and  be  respectively,  that  belong  to  C.  Since  C  is  convex,  the  line 
segment  />«»  Ae  is  contained  ia  C.  Also,  the  line  segment  P4t  P*e  Is  Inside  the  triangle  abc  and  must  intersect  the  diagonal 
bd  at  a  point  between  6  and  Pcmtrtid  (which  is  on  the  diagonal  ac) ,  say  Pi.  Thus  there  is  one  point  on  the  diagonal  bd 
between  6  and  Pcentroid  that  is  also  In  C.  Similarly  we  can  show  that  there  is  one  point  on  the  diagonal  bd  between  d  and 
Pccntroid,  say  Pi,  that  is  also  in  C.  Therefore,  since  C  is  convex,  P„„,r«a  is  part  of  C. 

End  of  Proof. 

This  theorem  suggests  that  the  neighborhood  of  the  centroid  is  a  good  place  to  look  for  possible  integer  solutions.  We 
use  this  fact  in  the  algorithm  presented  later  in  this  section.  Furthermore,  if  none  of  the  closest  integer  point  neighbors  of 
the  centroid  are  in  the  convex  region,  the  convex  region  can  only  be  of  a  special  shape.  Specifically  we  have  the  following 
results: 

Lemma  5  Let  R  (vertices  r  i  rirju)  and  Rf  (vertices  P2Pj  r^l  be  rectangles  with  parallel  sides  such  that  Rf  is  enclosed 
inside  R.  Let  p  be  a  point  inside  Rf.  A  set  of  line  segments  that  connect  a  point  on  each  side  of  R  top  must  intersect  at 
least  two  sides  of  Rf. 

Proof: 

We  will  use  Figure  3  for  illustration.  Any  line  segment  connecting  p  to  a  point  on  a  side  of  R  must  intersect  at  least 
one  side  of  the  rectangle  Rf  since  all  of  Pis  entirely  outside  Rf  and  pis  inside/?.  Suppose  r^rj  is  the  only  side  of  Rf  that 
the  set  of  line  segments  connecting  p  to  each  side  of  R  intersect  Now  p  is  a  point  between  parallel  lines  passing  through 
pjP2  and  P4P3  and  a  line  segment  connecting  a  point  on  rirj  to  p  cannot  intersect  the  line  through  P4P}.  Thus  r'4r '3  cannot 
be  the  only  side  of  Rf  that  the  set  of  line  segments  connecting  p  to  each  side  of  R  Intersect.  Hence,  there  must  be  at  least 
two  such  sides,  and  that  proves  the  lemur- 

End  of  Proof. 

Theorem  6  Let  C  be  a  convex  region  in  two  dimensions.  Let  R  be  the  smallest  rectangle  with  sides  parallel  to  the  axes 
that  encloses  C.  Assume  that  the  centroid  of  R  is  not  an  integer  point  and  let  Rf  be  the  smallest  rectangle! square)  with 


it 
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Figure  3:  Illustration  for  Lemma  5 


integer  vertices  that  includes  the  centroid  ofFLlf  none  of  vertices  of  R'  belong  to  C,  and  R  is  more  than  4  integer  units 
in  length  and  width,  then  C  intersects  exactly  two  sides  of  Rf. 

Proof: 

Let  R  be  the  rectangle  rlrjrjr4  ia  Figure 4(i)  that  tightly  bounds  a  convex  region  C  (not  shown  in  the  figure).  R'  is 
the  rectangle  shown  as  and  the  centroid  of  R  is'  an  interior  point  of  R! .  The  sides  of  Rf  are  extrapolated  until 

they  meet  a  side  of  ft. 

Since  R  tighly  bounds  C,  there  is  at  least  one  point  belonging  to  C  on  each  side  of  rectangle  R.  Also  there  is  at  least 
one  point  belonging  to  C  inside  Rf  (centroid  of  R).  Since  C  is  convex,  there  is  a  line  segment  from  a  point  on  each  side 
of  A  to  a  point  inside  Rf  which  is  completely  contained  in  R.  From  Lemma  5,  C  must  intersect  at  least  two  sides  of  Rf. 

We  now  prove  that  C  can  intersect  at  most  two  sides  of  Rf.  Suppose  C  intersects  three  sides  of  Rf.  Without  loss 
of  generality  we  assume  the  sides  are  r£r£,  and  r1^.  We  represent  this  by  drawing  these  sides  with  thick  lines  in 
Figure  4(ii).  In  Figure  4  we  will  convert  every  line  segment  that  we  prove  has  a  point  belonging  to  C,  to  thick  lines. 
Every  line  segment  which  is  proved  not  to  contain  any  points  in  C  is  marked  by  dashed  lines. 

We  have  assumed  that  there  is  at  least  one  point  on  the  line  segment  r^r/,  that  belongs  to  C.  If  any  point  on  the  line 
segment  rja  belongs  to  C,  then  the  point  must  belong  to  C,  since  C  is  convex.  But  r\  is  a  vertex  of  Rf  and  does  not 
belong  to  C.  Thus  we  conclude  that  the  line  segment  r{ a  does  not  intersect  C.  Similarly  we  can  show  that  line  segments 
rja,  rj/,  rji,  rje,  rjj  and  r^d  cannot  Intersect  C.  These  inferences  are  shown  in  Figure  4(iii). 

If  there  is  any  point  of  C  inside  the  rectangle  gr'Jr4,  then  a  line  segment  from  that  point  to  a  point  inside  Rf  must  be 
contained  in  C.  But  such  a  line  segment  must  intersect  gRA  or  r'/,  and  we  have  established  that  these  line  segments  do 
not  intersect  C.  Thus  there  can  be  no  points  inside  the  rectangle  gRJu  that  belong  to  C.  In  particular,  line  segments  /r4 
and  ug  do  not  intersea  C.  Similarly  we  can  show  that  line  segments  dr3  and  rje  do  not  contain  any  points  belonging  to 
<7.  This  is  shown  in  Figure  4(iv). 

Since  A  tightly  bounds  C,  there  must  be  at  least  one  point  on  nr4  that  is  in  C.  Since  we  have  proved  that  1  ine  segments 
r3e  and  /r4  do  not  contain  any  points  belonging  to  C,  the  line  segment  e /  must  intersea  C,  as  shown  in  Figure  4(v). 

Now  since  the  sides  of  Rate  more  than  4  integer  units  (or  4  times  the  sides  of  Rf),  line  segments  prj,  hr{,  r'J,  r'e 
are  all  longer  than  one  integer  unit  (or  longer  than  sides  of  Rf).  Also  sides  of  Rf  and  line  segments  ef  and  gh  are  of  unit 
length. 

It  can  be  seen  that  any  line  through  a  point  on  the  line  segment  ef  that  does  not  intersea  the  line  segment  fr'4  cannot 
intersect  the  tine  through  points  h  and  c  at  a  point  more  than  1  unit  left  (towards  h)  of  point  ,  and  hence  cannot  intersect 
the  tine  segment  gh  (see  pq  in  Figure  4(v».  Thus  there  cannot  be  points  belonging  to  C  on  both  line  segments  ef  and 
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Figure  4:  Intersection  of  a  convex  region  and  its  enclosing  rectangle 


gh,  since  a  line  segment  connecting  those  points  would  have  to  intersect  the  line  segment  fr^  and  we  have  proved  that 
the  line  segment  /rj  does  not  intersect  C.  Since  we  have  also  proved  that  the  line  segment  ef  has  at  least  one  point 
belonging  to  C,  the  line  segment  gh  cannot  intersect  C.  Similarly  we  can  show  that  the  line  segment  cd  cannot  intersect 
C.  Thus  we  reach  the  situation  shown  in  Figure  4{v). 

Now  since  the  side  r4ri  of  rectangle  R  must  have  at  least  one  point  that  is  in  C  and  the  line  segment  r4h  has  been 
shown  to  not  have  any  points  belonging  to  C,  the  line  segment  rt  h  must  intersect  C.  Similarly  we  can  show  that  the  line 


segment  rjc  must  intersect  C  (see  Figure  4(vi)).  Thus  there  must  be  a  line  segment  from  a  point  on  the  line  segment  n  h 
to  a  point  on  the  Use  segment  r3c  which  is  contained  in  C.  But  every  such  line  segment  must  intersect  the  line  segment 
ar{  (and  br'j)  which  we  have  proved  does  not  intersect  C.  Therefore  we  have  a  contradiction. 

Hence  the  convex  region  C  cannot  intersect  three  or  more  sides  of  Rf.  Since  we  have  already  shown  that  it  must 
intersect  at  least  two  sides  of  Rf.  the  convex  region  C  must  intersect  exactly  two  sides  of  R'. 

End  of  Proof. 

Furthermore,  we  find  that  a  convex  region  that  satisfies  the  requirements  of  this  theorem  must  be  located  around  one 
of  the  diagonals  of  the  bounding  rectangle.  Specifically  we  have  the  following  theorem: 

Theorem  7  Let  C,  R  and  Rf  be  defined  as  in  Theorem  6.  Let  R\,  R%  R)  and  R  be  the  rectangles  formed  by  the 
extensions  of  the  sides  of  Rf  and  the  rectangle  R  as  shown  in  Figure  5.  Then  under  the  conditions  of  Theorem  6.  exactly 
two  diagonally  opposite  rectangles  R  intersect  with  C. 


U  ft  rj 


Figure  5:  Illustration  for  Theorem  7 


Proof: 

We  first  show  that  at  least  two  of  the  rectangles  R  intersect  C.  Suppose  none  of  the  rectangles  R  intersect  C.  Since 
R  strictly  bounds  C,  all  sides  of  R  intersect  C.  Thus  C  must  intersect  each  of  the  segments  ab,  cd,  ef  and  gh,  since  those 
are  the  onty  parts  of  the  sides  of  R  that  do  not  belong  to  any  of  Ri.  But  a  line  segment  connecting  any  point  on  a6  to  any 
point  on  be  must  intersect  Ri,  which  contradicts  that  none  of  the  rectangles  R  intersects  C.  Hence  there  must  be  at  least 
one  rectangle  Ri  that  intersects  C. 

Suppose  R{  is  the  only  rectangle  R  that  intersects  C.  Then  C  must  intersect  the  line  segments  cd  and  ef  since  tnose 
are  the  only  parts  of  the  sides  rjr3  and  r3r«  respectively  that  do  not  belong  to  any  of  the  R  assumed  not  to  intersect  C. 
But  a  line  segment  connecting  any  point  on  cd  to  any  point  on  ef  must  intersect  the  rectangle  Rj,  which  contradicts  that 
i?i  is  the  only  rectangle  that  intersects  C.  We  conclude  that  at  least  two  of  the  rectangles  R  must  intersect  C. 

We  now  show  that  not  more  than  two  of  the  rectangles  R  can  intersect  C.  Suppose  at  least  ihree  of  the  rectangles  R 
intersect  C.  Let  the  rectangles  be  Rt,  Ri  and  Rj.  Let  pt,  pi  and  be  points  belonging  to  C  in  the  rectangles  Rt,  Ri 
and  R)  respectively.  It  is  easy  to  see  from  Figure  5  that  point  r\  is  inside  the  triangle p\pips,  and  hence  inside  C  which 
is  a  contradiction.  Hence  we  can  have  at  most  two  rectangles  R  that  intersect  C.  We  conclude  that  exactly  two  of  the 
rectangles  R  intersect  C. 
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Figure  6:  Illustration  for  Theorem  7 


We  now  show  that  the  two  rectangles  ft,  that  intersect  C,  must  be  diagonally  opposite.  Suppose  two  adjacent 
rectangles  fti  and  Ri  intersea  C.  Then  ft3  and  ft*  cannot  intersea  C,  since  exactly  two  of  ft,  intersect  C.  The  line 
segment  ef  must  intersea  C  since  that  is  the  only  part  of  the  side  nr 4  that  is  not  in  the  rectangles  ft3  or  ft,.  Let  p,  and 
P2  be  points  in  the  rectangles  Ri  and  R2  respectively  that  belong  to  C,  as  shown  in  the  Figure  6.  A  line  connecting  pt  and 
pi  will  intersea  the  segment  br^.  La  the  point  of  intersection  be  point  q\.  Also,  since  C  is  convex  and  tightly  bounded 
by  ft,  there  must  be  a  line  segment  contained  in  C  that  connects  a  point  on  the  side  rjr»  to  a  point  on  the  line  segment  /e. 
Since  the  line  segment  ab  is  shorter  than  the  line  segment  br2,  and  the  line  segment  er'2  is  longer  than  the  line  segment 
rjb,  any  line  segment  connecting  a  point  on  the  side  nn  to  a  point  on  the  line  segment  ef  must  intersect  line  segment 
rje,  say  at  point Since  qi  and  q2  belong  to  C  and  lies  on  the  line  connecting  them,  rj  must  belong  to  C.  This  is  a 
contradiction  to  the  conditions  of  this  theorem.  Hence  two  adjacent  Ri  cannot  intersea  C. 

Hence,  exactly  two  diagonally  opposite  Ri  must  intersea  C. 

End  of  Proof. 

We  conclude  that  if  a  convex  region  (larger  than  a  certain  size)  contains  an  integer  point,  then  either  an  integer  point 
close  to  the  centroid  of  its  hounding  rectangle  is  inside  the  convex  region,  or  the  convex  region  region  can  be  divided  into 
two  smaller  regions,  at  least  one  of  which  must  contain  an  integer  poinL  Based  on  these  results,  we  have  the  following 
procedure  to  determine  if  a  given  convex  region  contains  any  integer  points, 
procedure:  FindIntegers(C,  R ) 

input:  Convex  region  C  defined  by  inequalities  in  two  variables,  say  z  and  y.  Rectangle  R  tightly  bounding  C. 
output:  Boolean  value  representing  whether  an  integer  solution  exists  in  the  region. 

1.  If  R  covers  more  than  4  integers  in  both  dimensions,  GOTO  step  2.  Otherwise,  without  loss  of  generality,  say  r 
is  the  variable  which  can  take  at  most  4  integer  values.  For  each  possible  integer  value  of  z.  determine  the  range 
of  values  y  can  take.  If  an  integer  is  found  in  the  value  range  of  y,  return  TRUE,  else  return  FALSE.  This  is  an 
exhaustive  search  in  a  small  space. 

2.  Find  the  centroid  of  R.  If  it  is  an  integer  point,  return  TRUE.  La  Rf  be  the  smallest  square  with  integer  points  as 
corners  enclosing  the  centroid  of  R.  If  any  of  the  four  corner  points  of  R!  is  part  of  C,  return  TRUE.  Otherwise 
GOTO  step  3. 
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C :  Convex  region  under  consideration. 

R :  Smallest  rectangle  enclosing  C 

Rf :  Smallest  integer  square  enclosing  centroid  of  R 

R\,Ri\  Bounding  rectangles  for  the  new  problem 

Figure  7:  Nature  of  Convex  Regions  with  no  points  in  the  vicinity  of  the  Centroid  of  the  Bounding  Rectangle 


3.  Find  the  two  sides  of  Rf  that  intersect  C.  Determine  the  two  rectangles  Ri  and  as  shown  in  Figure  7.  Define 
Ci  and  Ci  by  adding  the  new  boundaries  of  Ri  and  Rj  respectively,  to  C  as  additional  inequality  constraints.  Find 
the  new  strictly  bounding  rectangles  R\  and  for  C\  and  Ci  respectively  using  Fourier-Motzldn  elimination. 
Return  (FmdIntegers(C|,  /?,)  OR  FindlntegersCCj,  R^)).  ■ 

complexity:  This  procedure  divides  the  problem  into  two  smaller  problems  that  are  half  the  size.  The  division  stops  at 
problem  size  4.  Thus  the  maximum  number  of  times  the  above  steps  may  have  to  be  repeated  is  n/4,  where  n  is  the 
length  of  the  bounding  rectangle  for  the  original  convex  region  C. 

end  of  procedure 

4.2  Solution  for  n  Dimensions 

The  results  obtained  for  2  dimensions  cannot  be  directly  extended  to  n  dimensions.  However,  intuition  suggests  that  the 
centroid  of  the  bounding  rectangle  is  still  a  good  position  to  potentially  establish  an  integer  point  We  use  this  to  divide 
and  shrink  the  region  for  analysis,  until  a  solution  is  found  or  the  problem  is  reduced  to  2-dimensions, 
procedure;  NdImFlndIntegers(n,C,  ft) 

input:  A  set  of  inequalities  with  n  variables  representing  a  closed  convex  region  C.  Values  A,”"  and  X’nin  tori  =  1  to  n 
denoting  an  n-rectangle  R  strictly  enclosing  C. 

output :  Boolean  value  representing  whether  an  integer  solution  exists  in  the  region. 


1 .  If  range  in  any  dimension  does  not  contain  an  integer,  no  integer  solutions  can  exist  so  terminate  returning  FALSE. 
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If  the  product  of  the  the  ranges  is  a  small  Integer,  perform  an  exhaustive  search  for  solutions  in  the  integer  space. 
If  the  number  of  dimensions  is  two.  Return  (FindIntegers(C,  R)).  Otherwise  go  to  step  2. 

2.  Determine  the  centroid  of  the  given  n-rectangle  R.  Let  Rf  be  the  smallest  ^-rectangle  with  all  integer  vertices  that 
contains  the  centroid  of  R.  If  any  of  the  integer  vertices  of  K  belong  to  the  convex  region  C,  return  TRUE,  else 
GOTO  step  3. 

3.  Let  X{  be  the  dimension  in  which  R  contains  the  smallest  number  of  integers.  If  the  number  of  integers  is  greater  than 
two,  GOTO  step  4.  Otherwise  obtain  (at  most)  two  new  convex  regions  in  n  -  1  dimensions  by  giving  occurrences 
of  Xi  fixed  integer  values  within  this  range.  Suppose  Cj  and  Cj  are  the  new  regions  obtained.  Find  the  smallest 
enclosing  n-rectangles  Rj  and  Hj.  Return  ( NdimFindIntegers(n- 1 ,  Cj ,  Ri )  OR NdimFmdIntegers(n- 1 ,  C2 ,  R2)). 

4.  Let  *i  and  xi  be  the  two  values  in  the  center  of  the  range  of  values  with  z\  <  z%  that  Xi  can  take.  Add  additional 
constraints  X,-  <=  zt  and  Xi  >=  zj  to  C  to  get  two  new  convex  regions  Cj  and  C2  respectively.  Find  their  smallest 
bounding  n-rectangles  Rj  and  R2.  Return  ( NdimFindlntegersfn.Cj,  R{)  OR  NdimFindlntegersfn.Cj,  R2)). 

end  of  procedure 

The  steps  of  these  solution  procedures  are  organized  so  that  most  common  cases  can  be  handled  quickly,  but  the  less 
common  cases  are  also  handled.  It  is  not  possible  to  do  an  accurate  complexity  analysis  of  this  solution  procedure  since  it 
is  dependent  on  numerical  values.  However,  we  expect  the  procedure  to  be  efficient  for  problems  with  very  few  variables 
and  small  ranges,  which  is  what  we  expect  in  our  applications. 

43  Handling  Infinite  Regions 

The  solution  procedures  so  far  assumed  that  the  ranges  for  all  variables  are  finite.  This  may  not  be  true  in  general.  An 
infinite  convex  region  would  include  zero  or  infinite  integer  points.  It  is  obvious  that  an  infinite  convex  region  bounded 
by  diverging  byperplanes  would  include  infinite  integer  points.  A  necessary  condition  for  an  infinite  convex  region  to 
have  no  integer  points  is  that  there  are  parallel  hyperplanes  bounding  the  region  (with  an  infinite  thin  strip  with  no  integer 
points  between  them).  We  have  the  following  lemma; 

Lemma  8  An  infinite  size  convex  region  must  have  infinite  integer  points  unless  at  least  two  of  the  hyperplanes  defining 
it  are  parallel 

For  our  analysis,  we  check  if  any  of  the  hyperplanes  defining  the  convex  region  are  parallel.  If  no  two  hyperplanes 
are  parallel,  no  further  analysis  is  required.  If  parallel  hyperplanes  do  exist,  we  generate  heuristic  values  for  bounds 
whenever  actual  values  are  not  available.  A  heuristic  value  that  we  use  for  bounds  is  the  largest  coefficient  in  the  system 
of  inequalities.  This  is  supported  by  the  fact  that  the  occurrence  of  integers  inside  infinite  convex  regions  is  usually 
cyclic  along  the  axes  in  which  it  is  infinite,  and  the  lengths  of  cycles  are  determined  by  the  magnitude  of  the  integers  in 
the  inequality  set.  A  more  exact  analysis  could  be  performed  to  determine  a  better  defined  bound,  but  we  do  not  have 
evidence  that  would  suggest  that  our  heuristic  is  inadequate. 

5  Related  Work 

Exact  data  dependence  analysis  has  been  addressed  by  several  research  groups  in  the  recent  past.  We  compare  our  method 
to  the  work  of  Maydan  et.  al.  [6J,  Omega  test  [7]  and  Power  test  (11],  The  first  step  in  our  solution  procedure,  that  is 
transformation  of  equality  constraints  to  inequality  constraints  (which  may  yield  a  solution  in  some  cases),  is  also  used 
in  the  approaches  stated  above.  The  difference  in  the  methodologies  used  is  discussed  in  section  2.4.  The  second  step 
that  we  apply  is  Fourier-Motzkin  elimination  to  determine  the  existence  of  a  real  solution  to  the  transformed  system.  All 
the  other  methods  also  use  Fourier-Motzkin  elimination  when  needed.  Maydan  et.  al.  examine  the  form  of  the  equation 
system  at  this  stage  and  check  if  any  of  a  set  of  simple  special  tests  is  applicable,  before  employing  Fourier-Motzkin 
elimination.  In  power  test,  a  small  approximation  is  introduced  to  enable  an  efficient  solution  methodology,  and  hence 
the  results  may  be  inexact.  In  Omega  test,  a  modified  version  of  Fourier-Motzkin  elimination  is  used,  which  can  prove 
the  existence  of  an  integer  solution  in  many  situations.  We  have  not  specifically  addressed  the  best  way  of  performing 
this  step,  but  many  optimizations  are  certainly  possible.. 
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The  unique  feature  of  our  approach  is  the  Anal  step,  which  is  reached  if  none  of  the  above  two  steps  are  able  to  prove 
or  disprove  the  existence  of  an  integer  solution.  Maydan  et.  al.  suggest  using  a  branch  and  bound  strategy  at  this  stage. 
In  Omega  test,  the  final  equation  system  is  solved  by  successively  eliminating  variables  by  building  a  new  system  for 
every  feasible  value  of  the  variable.  We  employ  a  sophisticated  search  strategy  based  on  geometric  properties  of  convex 
regions.  In  particular,  the  fact  that  an  integer  solution  in  a  closed  convex  region  is  likely  to  be  close  to  the  cenuoid  of  the 
region,  is  used  to  guide  the  search  procedure.  If  a  step  does  not  yield  a  solution,  it  generates  a  pair  of  two  smaller,  more 
constrained  convex  regions  as  the  possible  solution  spaces. 


6  Conclusions 

We  have  presented  a  methodology  for  determining  whether  a  non-negative  integer  solution  for  a  system  of  equations 
exists.  Our  experience  suggests  that  the  method  is  suitable  and  sufficient  for  equation  systems  that  need  to  be  solved 
in  a  compiler  for  parallel  languages.  The  research  was  motivated  by  the  problem  of  determining  if  synchronization  in 
a  parallel  program  is  sufficient  to  ensure  sequentially  consistent  results  [9].  It  is  particularly  important  to  have  an  exact 
solution  when  examining  race  conditions  in  parallel  programs,  since  the  existence  of  an  exact  solution  is  used  to  guarantee 
deterministic  parallel  execution. 

However,  the  approach  is  general  and  can  be  used  for  data  dependence  analysis  and  other  small  problems  with  similar 
characteristics.  Although  the  problems  involved  are  NP-hard,  in  practice  we  can  solve  almost  all  real  instances  without 
much  trouble,  because  of  practical  limits  on  their  sizes. 
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